NEFSC Exploratory Data Analysis

Exporatory Data Analysis of NEFSC Bottom Trawl Survey Data

Exploration of spatial and temporal patterns in abundance, and bodymass of fishes from the Northeast groundfish survey.

Build code containing data wrangling and conversions can be accessed here.

Data loaded for this markdown is pulled directly from the {targets} pipeline for data consistency with a developing codebase.

####  Load NEFSC Groundfish Data  ####

# Loading from targets:
# Load the area-stratified biomass/abundances that used species cutoffs
withr::with_dir(rprojroot::find_root('_targets.R'), 
                tar_load(nefsc_stratified))


# Do some text formatting
nefsc_stratified <- nefsc_stratified %>% 
  mutate(
    id = as.character(id),
    season = ifelse(season %in% c("spring"), "Spring", "Fall"),
    season = factor(season, levels = c("Spring", "Fall"))
  )

# Run Summary Functions
ann_means <- ss_annual_summary(nefsc_stratified, include_epu = F)
seasonals <- ss_seasonal_summary(nefsc_stratified, include_epu = F) 

# bind them so you can facet
summs <- bind_rows(ann_means, seasonals) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))


#drop haddock to see if that changes the bump
haddock_ann  <- ss_annual_summary(filter(nefsc_stratified, comname != "haddock"))
haddock_seas <- ss_seasonal_summary(filter(nefsc_stratified, comname != "haddock"))
no_haddockn  <- bind_rows(haddock_ann, haddock_seas) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))

Temporal Trends

Patterns in annual abundance or CPUE of High Profile Species at annual and seasonal scales.

Total Biomass

Biomass of a given species from a station is a function of the adjusted number caught, their lengths, and the corresponding weights at-length using the species specific length-weight relationship coefficients. Sex-specific coefficients and seasonal growth coefficients are used when available.

Total biomass in the following plots represent the total of all individual biomass estimates from all species with available growth coefficients, for each year and within each survey season.

L-W Derived

Species specific length-weight relationships allow us to look at where biomass is across different lengths. When aggregating up to total biomass using these expected relationships we get this timeseries.

# Annual Timeline
bio_1 <- ann_means %>% 
  ggplot(aes(est_year, lw_biomass_kg, 
             color = season, 
             linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, lw_biomass_kg, 
                  color = season, 
                  linetype = "Haddock Removed")) +
    scale_y_continuous(labels = scales::comma_format()) +
    scale_color_manual(values = "black") +
    labs(x = "", 
         y = "Total Length-Weight Biomass (kg)",
         linetype = "",
         color = "Season")


# Seasonal Timeline
bio_2 <- seasonals %>% 
  ggplot(aes(est_year, lw_biomass_kg, color = season)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", 
         y = "Total Length-Weight Biomass (kg)", 
         color = "Season",
         caption = "Biomass totals from actual catch, not stratified to represent survey areas.")

# Stack them with patchwork
bio_1 / bio_2

Gross Shipboard Weight

As a comparison, we can also get the total biomass from the weights recorded on the ship. Not all fishes are weighed individually, but total biomass for a species is recorded for each station.

We would expect with accurate length-weight relationships that this plot should resemble the length-weight time line.

ship_1 <- ann_means %>% 
  ggplot(aes(est_year, fscs_biomass_kg, color = season, linetype = "All Species")) +
    geom_line(group = 1) +
    geom_line(data = haddock_ann, 
              aes(est_year, fscs_biomass_kg, color = season, linetype = "Haddock Removed")) +
    scale_color_manual(values = "black") +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Survey Biomass (kg)")

ship_2 <- seasonals %>% 
  ggplot(aes(est_year, fscs_biomass_kg, color = season)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Survey Biomass (kg)", caption = "Biomass from survdat$biomass.")

ship_1 / ship_2

Area Stratification

Area stratified totals for abundance and biomass are calculated according to the following formulas in accordance with code found here, courtesy of Sean Lucey

The stratum area ratio is the area of strata: \(1,2,3, ...s\) relative to the total area of all strata sampled that year (only stratum with data we use is included).

Alternatively, catch can be similarly weighted acccording to the effort and areas of the regional Ecological Production Units. Both of the following are done for each year to account for annual variation in coverage.

\[stratum~ratio = area(stratum_s)~ / total~area\]

or

\[epu~ratio = area(epu_s) / total~area\]

Next, using the amount of tows in a given year for a specific strata/epu, we get the amount of individuals across that group and divide it by the effort. This is then scaled by \(stratum~ratio\) to scale it by how abundance for that stratum reflects the total area of the surey.

For Strata: \(1,2,3, ...s\)
For years: \(1,2,3, ...i\)

\[stratum~weighted~catch = stratum~ratio_s * (\frac{~abundance_{is}~}{~n~stations_{is}}~),\]

The stratum weighted catch (abundance or biomass) is the average catch rate for a given area, weighted by how much the area that stratum is relative to the overall survey area that year.

The stratum weighted catch is then multiplied by the total survey area for that year, divided by the area of a standard tow by the RV Albatross (0.0384 square km). That value can be adjusted to reflect differences in catchability (q) if available, otherwise q is explicitly assumed to be 1. This gives the total abundance or biomass for a given size class across the entire strata based on its catch rate that year and the area that the stratum represents.

\[area~expanded~abundance =(stratum~weighted~catch * \frac{total~area}{.0384~square~km})~/~q\]

Then to get the expected biomass for a given size class using the length weight coefficients, the expanded total abundance for that size class just needs to be multiplied by what the mean bodymass would be for that length fish and its species-specific coefficients.

\[area~expanded~biomass = area~expanded~abundance~*~size~specific~biomass\]

Abundance

When expanding the catches out to the total area the survey represents we can get this time line for annual abundances.

abund_1 <- ann_means %>% 
  ggplot(aes(est_year, total_survey_abund)) +
  geom_line(color = "gray50", linetype = 3) +
  geom_line(aes(color = `Research Vessel`), show.legend = F) +
  scale_y_continuous(labels = scales::comma_format()) +
  scale_color_gmri() +
  scale_x_continuous(labels = scales::breaks_pretty()) +
  labs(x = "", y = "Total Fish from Survey") 

abund_2 <- ann_means %>% 
  ggplot(aes(est_year, strat_abundance_s/1000000)) +
  geom_line(color = "gray50", linetype = 3) +
  geom_line(aes(color = `Research Vessel`)) +
  scale_y_continuous(labels = scales::comma_format()) +
  scale_color_gmri() +
  scale_x_continuous(labels = scales::breaks_pretty()) +
  labs(x = "", y = "Area Extrapolated Abundance \n(millions)",
       caption = "Area stratification weights from nmfs trawl stratum, not ecological production units.") +
  theme(legend.position = "bottom")


abund_1 / abund_2

Size Specific L-W Biomass

Biomass below uses catch at length, and the expected biomass at length, to get size specific catch/tows across all years for each stratum.

These catch rates are then expanded from catch per tow, to expected total across the entire area for each strata using the difference in relative area among them and the expected area covered by a normal tow.

lw_strat_1 <-  ggplot(ann_means) +
    geom_line(aes(est_year, lw_strat_biomass_s/1000000, 
                  color = season)) +
   scale_color_manual(values = "black") +
    scale_y_continuous(labels = scales::comma_format()) +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    labs(x = "", 
         y = "Area Expanded Biomass\n l-W Regressions (million kg)",
         color = "") +
  theme(legend.position = "bottom")

lw_strat_2 <- ggplot(seasonals) +
    geom_line(aes(est_year, lw_strat_biomass_s/1000000, color = season)) +
    scale_y_continuous(labels = scales::comma_format()) +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    labs(x = "", 
         y = "Area Expanded Biomass\n L-W Regressions (million kg)", 
         color = "")+
  theme(legend.position = "bottom")

lw_strat_1 / lw_strat_2

Total Species Biomass

The totals below use the survdat$biomass column for size-agnostic biomass catch rates. Uses strata specific biomass catch rates which are then expanded out to be representative of total survey areas.

strat_1 <-  ggplot(ann_means) +
    geom_line(aes(est_year, fscs_strat_biomass_s/1000000, color = season)) +
    scale_color_manual(values = "black") +
    scale_y_continuous(labels = scales::comma_format()) +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    labs(x = "", 
         y = "Area Expanded Biomass \n survdat$BIOMASS (million kg)",
         color = "Season") +
  theme(legend.position = "bottom")

strat_2 <-  ggplot(seasonals) +
    geom_line(aes(est_year, fscs_strat_biomass_s/1000000, color = season)) +
    scale_y_continuous(labels = scales::comma_format()) +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    labs(x = "", 
         y = "Area Expanded Biomass \n survdat$BIOMASS (million kg)",
         color = "Season") +
  theme(legend.position = "bottom")

strat_1 / strat_2

CPUE

Biomass per station in the following figures represents the total biomass divided by the effort in number of stations.

There has been no correction for differences in physical stratum area in these figures, only station effort within them for a given year.

cpue_1 <- ann_means %>% 
  ggplot() +
    geom_line(aes(est_year, lw_biomass_per_station, 
                  color = season, linetype = "All Species")) +
    geom_line(data = haddock_ann, 
              aes(est_year, lw_biomass_per_station, 
                  color = season, linetype = "Haddock Removed")) +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    scale_color_manual(values = "black") +
    labs(x = "", 
         y = "LW Biomass / Station\n(kg)",
         linetype = "",
         color = "Season") +
  theme(legend.position = "bottom")

cpue_2<- seasonals %>% 
  ggplot(aes(est_year, lw_biomass_per_station, color = season)) +
    geom_line() +
    geom_vline(xintercept = 2008.5, linetype = 2, alpha = 0.5) +
    labs(x = "", 
         y = "LW Biomass / Station\n(kg)", 
         color = "Season") +
  theme(legend.position = "bottom")

cpue_1 / cpue_2

Species Guilds

Relative contribution to overall biomass by species group based on management and/or life history: Demersal, Elasmobranch, Groundfish, Pelagic

# Species Class Bins Data Prep
class_totals <- nefsc_stratified %>% 
  group_by(est_year, spec_class) %>% 
  summarise(`Total Class Biomass` = sum(sum_weight_kg),
            .groups = "keep") %>% 
  ungroup() %>% 
  left_join(ann_means, by = "est_year") %>% 
  mutate(`Percent of Annual Biomass` = (`Total Class Biomass` / lw_biomass_kg) * 100 ,
         `Species Guild` = spec_class,
         `Species Guild` = ifelse(is.na(`Species Guild`), "Unclassified", `Species Guild`),
         Year = as.character(est_year))

Proportion of Total Catch

Class Totals

# Class Proportions
class_totals %>% 
  ggplot(aes(x = fct_rev(Year), y = `Percent of Annual Biomass`, fill = `Species Guild`)) + 
    geom_col(color = "white", size = 0.2) +
    geom_vline(xintercept = "2008", linetype = 2, color = "gray25") +
    labs(x = "", y = "Proportion of Total Annual Biomass") +
    coord_flip() +
    scale_fill_gmri() +
    scale_y_continuous(position = "right") +
    guides("fill" = guide_legend(title.position = "top", title.hjust = 0.5)) +
    theme(legend.position = "bottom")

Plotting Individual Guilds

guild_plots <- nefsc_stratified %>% 
  split(.$spec_class) %>% 
  imap(function(guild, guild_name){
    
    
    year_df <- guild %>% 
      split(.$est_year) %>% 
      imap_dfr(function(year_group, year){
        #quantile breaks
        probs <- c(0.5, 0.25, 0.5, 0.75, 0.95 )
        
        # length quantiles
        length_quantiles <- Hmisc::wtd.quantile(
            x = year_group$length, 
            weights = year_group$numlen_adj, 
            probs = ) %>% 
          as.data.frame() %>% 
          rownames_to_column() %>% 
          setNames(c("Quantile", "length")) %>% 
          mutate(Quantile = factor(probs),
                 year = year)
        
        # weight quantiles
        weight_quantiles <- Hmisc::wtd.quantile(
            x = year_group$ind_weight_kg, 
            weights = year_group$numlen_adj, 
            probs = probs) %>% 
          as.data.frame() %>% 
          rownames_to_column() %>% 
          setNames(c("Quantile", "weight (kg)")) %>% 
          mutate(Quantile = factor(probs),
                 year = year)
        
        guild_data <- left_join(length_quantiles, weight_quantiles, 
                               by = c("year", "Quantile")) %>% 
          mutate(year = fct_rev(year),
                 `Research Vessel` = ifelse(as.numeric(as.character(year)) > 2008, "HB", "AL"))
        
        return(guild_data)  
          
    
      })
    
    
    
    
    # Make key for better labels
    guild_key <- c("dem" = "Demersals", 
                   "el" = "Elasmobranchs", 
                   "gf" = "Groundfish", 
                   "pel" = "Pelagics")
    
    
    # build plot
    p1 <- ggplot(year_df, aes(y = length, x = factor(year))) +
      geom_rect(ymin = -Inf, 
                ymax = Inf, 
                xmin = "2008", 
                xmax = "2019", 
                fill = "gray90") +
      geom_line(aes(group = year), color = "gray50", alpha = 0.5) +
      geom_point(aes(color = Quantile)) +
      scale_color_gmri() +
      labs(subtitle = guild_key[guild_name], y = "Individual Length (cm)", x = NULL) +
      theme(legend.position = "none",
            axis.text.x = element_blank(),
            axis.ticks.x = element_blank())
    
    # build weight plot
    p2 <- ggplot(year_df, aes(y = `weight (kg)`, x = factor(year))) +
      geom_rect(ymin = -Inf, 
                ymax = Inf, 
                xmin = "2008", 
                xmax = "2019", 
                fill = "gray90") +
      geom_line(aes(group = year), color = "gray50", alpha = 0.5) +
      geom_point(aes(color = Quantile)) +
      scale_color_gmri() +
      labs(y = "Individual Weight (kg)", x = NULL) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) 
      
    patch_plot <- p1 / p2
    return(patch_plot)
    
    
  })

Pelagics

guild_plots$pel

Demersals

guild_plots$dem

Elasmobranchs

guild_plots$el

Groundfish

guild_plots$gf

Mean Individual Biomass

Check these across years and also by season

summs %>% 
  ggplot(aes(est_year, mean_ind_bodymass_lw * 1000, color = season)) +
    geom_line(show.legend = FALSE) +
    scale_color_manual(values = c(
      "Spring + Fall" = "black",
      "Spring" = as.character(gmri_cols("gmri blue")),
      "Fall"   = as.character(gmri_cols("orange")))) +
  labs(x = "", y = "Mean Individual Biomass (g)") +
  facet_grid(season~area)

Mean Individual Length

summs %>% 
  ggplot(aes(est_year, mean_ind_length, color = season)) +
    geom_line(show.legend = FALSE) +
    scale_color_manual(values = c(
      "Spring + Fall" = "black",
      "Spring" = as.character(gmri_cols("gmri blue")),
      "Fall"   = as.character(gmri_cols("orange")))) +
  labs(x = "", y = "Mean Individual Length (cm)") +
  facet_grid(season~area)

Top 5 Species by Weight

Since 2008 there has been a large increase in the amount of total biomass caught and the adjusted CPUE. As a way to look into the top species contributing to this I have pulled out the top 5 species for that year by total biomass caught as derived from their length weight relationships.

# Pull top 5 species from each year, by biomass
top_5 <- nefsc_stratified %>% 
  filter(est_year >= 2004) %>% 
  group_by(est_year, comname) %>% 
  summarise(n_fish = sum(numlen_adj),
            lw_biomass = sum(sum_weight_kg)) %>% 
  ungroup() %>% 
  split(.$est_year) %>% 
  map_dfr(function(yr_group){
    yr_group  %>% 
      slice_max(order_by = lw_biomass, n = 5) })
  

# Plot each year, re-ordered by biomass
t5_plots <- top_5 %>% 
  split(.$est_year) %>% 
  map(function(x){
    x %>% mutate(comname = fct_drop(comname),
                 comname = fct_reorder(comname, .fun = max, lw_biomass)) %>% 
    ggplot(aes(x = lw_biomass, y = comname)) +
      geom_point() +
      geom_segment(aes(yend = comname, x = 0, xend = lw_biomass)) +
      facet_wrap(~est_year) +
      scale_x_continuous(labels = scales::comma_format()) +
      labs(x = "Total Biomass (kg) - From LW Relationship", y = "")
})



# Assemble 4-panel plots displaying years in groups of 4
t1 <- ((t5_plots$`2004` | t5_plots$`2005`) / (t5_plots$`2006` | t5_plots$`2007`))
t2 <- ((t5_plots$`2008` | t5_plots$`2009`) / (t5_plots$`2010` | t5_plots$`2011`))
t3 <- ((t5_plots$`2012` | t5_plots$`2013`) / (t5_plots$`2014` | t5_plots$`2015`))
t4 <- ((t5_plots$`2016` | t5_plots$`2017`) / (t5_plots$`2018` | t5_plots$`2019`))

2004 - 2007

t1

2008 - 2011

t2 

2012 - 2015

t3

2016 - 2019

t4

Small Fish

Bodymass size cutoffs are used to plot the timelines of total fish biomass below the following thresholds for individual biomass: .0625, .125, .25, .5, 1kg. This is going to fail for the shipboard weights because that weight is divided equally across different measurements of fishes regardless of size.

.0625 kg

max_weight_plot <- function(size_cutoff){
  nefsc_stratified %>% 
    filter(ind_weight_kg < size_cutoff) %>% 
    group_by(est_year) %>% 
    summarise(total_biomass_lw = sum(sum_weight_kg),
              `Maximum Individual Biomass` = str_c("Individual Mass < ", size_cutoff, " kg"),
              .groups = "keep") %>%
    ggplot(aes(est_year, total_biomass_lw, color = `Maximum Individual Biomass`)) +
      geom_line() +
      scale_color_gmri() +
      scale_y_continuous(labels = scales::comma_format()) +
      labs(x = NULL, y = "Total Annual (kg)") +
    theme(legend.position = c(0.3, 0.9))
}


# .5 kg
max_weight_plot(.0625)

.125 kg

max_weight_plot(.125)

.25 kg

max_weight_plot(.25)

.5

max_weight_plot(.5)

1 kg

max_weight_plot(1)

Vessel Changes

Maximum size plots show a huge increase in biomass of small fishes corresponding with the vessel change from the RV Albatross to the RV Henry Bigelow.

To Take a look at what species contributed the most to this, I will maintain the maximum body size filter of .125 kg, and compare the difference in catch/tow for species across research vessel. The following plots will highlight the top & bottom ten species experiencing positive/negative changes to their catch rates. Temporal differences are ignored here, and may play a bigger role than vessel changes.

# annual biomasses for top 10 species, ordered by total biomass that year
# Minimum individual bodymass of 5k

# Compare total biomass and cpue across vessels
svp <- nefsc_stratified %>% 
    filter(ind_weight_kg < .125) %>% 
    group_by(svvessel, comname) %>% 
    summarise(bio = sum(sum_weight_kg),
              effort = n_distinct(id),
              cpue = bio / effort,
              .groups = "keep") %>% 
  pivot_wider(names_from = svvessel, 
              names_sep = "_",
              values_from = c(bio, effort, cpue)) %>% 
  mutate(cpue_change = (cpue_HB - cpue_AL) / cpue_AL) %>% 
  arrange(desc(cpue_change)) %>% 
  ungroup

Highlighting Positive Shifts in CPUE

# Highlight Positive Changes in CPUE
pos_change_plot <- svp %>% 
  slice_max(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(.f = comname, .x = cpue_change, .fun = max)) %>% 
  
  ggplot(aes(y = comname, x = cpue_change)) +
  geom_segment(aes(x = 0, xend = cpue_change, y = comname, yend = comname)) +
  geom_point() +
  labs(x = "Relative CPUE Change \n AL -> HB", y = "", subtitle = "Top Positive Shifts in CPUE")

# individual cpues
al_p <- svp %>% 
  slice_max(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(comname, cpue_AL, max)) %>% 
  ggplot(aes(y = comname, x = cpue_AL)) +
  geom_segment(aes(x = 0, xend = cpue_AL, y = comname, yend = comname)) +
  geom_point() +
  labs(x = "Biomass (kg) / Station", y = "", subtitle = "Albatross CPUE")

hb_p <- svp %>% 
  slice_max(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(comname, cpue_HB, max)) %>% 
  ggplot(aes(y = comname, x = cpue_HB)) +
  geom_segment(aes(x = 0, xend = cpue_HB, y = comname, yend = comname)) +
  geom_point() +
  labs(x = "Biomass (kg) / Station", y = "", subtitle = "Bigelow CPUE")


pos_change_plot | (al_p / hb_p)

Highlighting Negative Shift in CPUE

# Highlight Negative Shift in CPUE
neg_change_plot <- svp %>% 
  slice_min(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(comname, cpue_change, max)) %>% 
  ggplot(aes(y = comname, x = cpue_change)) +
  geom_segment(aes(x = 0, xend = cpue_change, y = comname, yend = comname)) +
  geom_point() + 
  scale_x_reverse() +
  labs(x = "Relative CPUE Change \n AL -> HB", y = "", subtitle = "Top Negative Shifts in CPUE")

# individual cpues
al_neg <- svp %>% 
  slice_min(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(comname, cpue_AL, min)) %>% 
  ggplot(aes(y = comname, x = cpue_AL)) +
  geom_segment(aes(x = 0, xend = cpue_AL, y = comname, yend = comname)) +
  geom_point() +
  labs(x = "Biomass (kg) / Station", y = "", subtitle = "Albatross CPUE")

hb_neg <- svp %>% 
  slice_min(n = 10, order_by = cpue_change) %>% 
  mutate(comname = fct_reorder(comname, cpue_HB, min)) %>% 
  ggplot(aes(y = comname, x = cpue_HB)) +
  geom_segment(aes(x = 0, xend = cpue_HB, y = comname, yend = comname)) +
  geom_point() +
  labs(x = "Biomass (kg) / Station", y = "", subtitle = "Bigelow CPUE")


neg_change_plot | (al_neg / hb_neg)

Big Fish

Same idea for the larger fish, but the opposite. A minimum cutoff for the individual bodymass. Below are the upper quantiles for individual body mass weighted by number at length.

# plotting function
min_weight_plot <- function(size_cutoff){
  nefsc_stratified %>% 
    filter(ind_weight_kg > size_cutoff) %>% 
    group_by(est_year) %>% 
    summarise(total_biomass_lw = sum(sum_weight_kg),
              `Minimum Individual Biomass` = str_c("Individual Mass > ", size_cutoff, " kg"),
              .groups = "keep") %>%
    ggplot(aes(est_year, total_biomass_lw, color = `Minimum Individual Biomass`)) +
      geom_line() +
      scale_color_gmri() +
      scale_y_continuous(labels = scales::comma_format()) +
      labs(x = NULL, y = "Total Annual (kg)") +
    theme(legend.position = c(0.3, 0.9))
}

1 kg

min_weight_plot(1)

2 kg

min_weight_plot(2)

5 kg

min_weight_plot(5)

10 kg

min_weight_plot(10)

Recent Upswings

Plots of minimum bodysize cutoff at 5kg and 10kg seem to have a big upswing in the last couple years, the following plot looks to highlight where the most biomass resides in those years with that same minimum bodymass cutoff of 5g.

# annual biomasses for top 10 species, ordered by total biomass that year
# Minimum individual bodymass of 5k
y_plots <- nefsc_stratified %>% 
    filter(ind_weight_kg > 5,
           est_year >= 2016) %>% 
    group_by(est_year, comname) %>% 
    summarise(total_biomass_lw = sum(sum_weight_kg),
              `Minimum Individual Biomass` = str_c("Individual Mass > ", 5, " kg")) %>%   
  split(.$est_year) %>% 
  map(function(x){
    x %>%  
      top_n(n = 20, total_biomass_lw) %>% 
      arrange(desc(total_biomass_lw)) %>% 
      mutate(comname = fct_reorder(comname, total_biomass_lw, mean))  %>% 
      ggplot(aes(total_biomass_lw, comname)) +
      geom_segment(aes(x = 0, xend = total_biomass_lw, comname, yend = comname)) +
      geom_point() +
      facet_wrap(~est_year) +
      labs(x = "Total LW Biomass (kg)", y = NULL)})

y_plots$`2016` | y_plots$`2017`

y_plots$`2018` | y_plots$`2019`

Cownose Rays

Big jumps in 2018 and 2019 seem to stem from catches of cownose rays. Further inspection shows that they boil down to two alarmingly high catches in the fall of each year. The table below details all instances of catch for cownose rays since 1990. Typically there have been 0-2 per year, with a deeply concerning increase in 2018-2019.

# How many cownose rays, and how suspicious is that number?
 cownose <- nefsc_stratified %>% 
   filter(comname == "cownose ray") %>% 
   group_by(comname, est_year, season, id) %>% 
   summarise(total_biomass_lw = sum(sum_weight_kg),
              num_individuals = sum(numlen_adj),
              n_occurrences   = n_distinct(id)) 

cownose %>% kable()

comname	est_year	season	id	total_biomass_lw	num_individuals	n_occurrences
cownose ray	1972	Fall	1972080531650	137.286135	16	1
cownose ray	1972	Fall	1972080551650	36.689235	4	1
cownose ray	1976	Fall	1976091121650	7.425116	1	1
cownose ray	2006	Fall	2006100691650	2.046616	1	1
cownose ray	2009	Fall	2009040721650	28.211572	2	1
cownose ray	2011	Fall	2011050531650	22.596350	1	1
cownose ray	2012	Fall	2012040861650	1.243419	1	1
cownose ray	2018	Fall	2018040781690	80.136882	45	1
cownose ray	2019	Fall	2019041521730	1771.707695	109	1

Size Spectra Bins

Using the same bodymass bins as in Kathy’s code cleaned up size spectra code.R we can take a look at how an overall size spectrum slope would look like for the two surveys. This will help us gain a coarse idea of where to put limits on the bodymass ranges based on the survey gear.

One thing that was brought up in my discussion with Andrew Edwards was that I should be standardizing the bin heights by their width because evenly spaced values on a log scale are increasingly large on the untransformed scale. For example for bins using log10() scale, the bin width from 2 < log10(x) < 3 is 900, and the next bin 3 < log10(x) < 4 has a width of 9,000.

# Set up bins for plotting normalized size specra

# COMMENTS
# bins increase in width as you work up on the log scale, 
# need to narmalize them by bin width so that they are comparable.

# NOTE: case_when is greedy so once something is assigned it is out of the pool of candidates
# In this case you want to count down    
nefsc_stratified <- nefsc_stratified %>% 
  mutate(
    lweight = log(ind_weight_kg),
    llen    = log(length),
    l10wt   = log10(ind_weight_kg))




 # Length bins for log() scale
nefsc_stratified <- nefsc_stratified %>% 
  mutate(
    # Bin Numbers
    lenbin2 = case_when(
      llen > 5    ~ 10,
      llen > 4.5  ~ 9,
      llen > 4    ~ 8,
      llen > 3.5  ~ 7,
      llen > 3    ~ 6,
      llen > 2.5  ~ 5,
      llen > 2    ~ 4,
      llen > 1.5  ~ 3,
      llen > 1    ~ 2,
      llen > 0.5  ~ 1,
      llen >   0  ~ 0),
    # Range of lengths for each bin
    len_range = case_when(
      lenbin2 == 10 ~ exp(5.5) - exp(5),
      lenbin2 == 9  ~ exp(5) - exp(4.5),
      lenbin2 == 8  ~ exp(4.5) - exp(4),
      lenbin2 == 7  ~ exp(4) - exp(3.5),
      lenbin2 == 6  ~ exp(3.5) - exp(3),
      lenbin2 == 5  ~ exp(3) - exp(2.5),
      lenbin2 == 4  ~ exp(2.5) - exp(2),
      lenbin2 == 3  ~ exp(2) - exp(1.5),
      lenbin2 == 2  ~ exp(1.5) - exp(1),
      lenbin2 == 1  ~ exp(1) - exp(.5),
      lenbin2 == 0  ~ exp(0.5) - exp(0)))




# Weight bins on log() scale
nefsc_stratified <- nefsc_stratified %>% 
  mutate(
    # Bin numbers
    wtbin = case_when(
      lweight > 8   ~ 10,
      lweight > 6   ~ 9,
      lweight > 4   ~ 8,
      lweight > 2   ~ 7,
      lweight > 0   ~ 6,
      lweight > -2  ~ 5,
      lweight > -4  ~ 4,
      lweight > -6  ~ 3,
      lweight > -8  ~ 2,
      lweight > -10  ~ 1,
      lweight < -10  ~ 0),
    # Weight Range for each bin
    wt_range = case_when(
      wtbin == 10 ~ exp(10) - exp(8),
      wtbin ==  9 ~ exp(8) - exp(6),
      wtbin ==  8 ~ exp(6) - exp(4),
      wtbin ==  7 ~ exp(4) - exp(2),
      wtbin ==  6 ~ exp(2) - exp(1),
      wtbin ==  5 ~ exp(1) - exp(.5),
      wtbin ==  4 ~ exp(.5) - exp(.25),
      wtbin ==  3 ~ exp(.25) - exp(.125),
      wtbin ==  2 ~ exp(.125) - exp(.0625),
      wtbin ==  1 ~ exp(.0625) - exp(.03125),
      wtbin ==  0 ~ exp(.03125) - exp(.015625)))




# Weight Bins on log10 scale
nefsc_stratified <- nefsc_stratified %>% 
  mutate(
    #Bin Numbers
    wt10bin =  case_when(
      l10wt > 3     ~ 21,
      l10wt > 2.5   ~ 20,
      l10wt > 2     ~ 19,
      l10wt > 1.5   ~ 18,
      l10wt > 1     ~ 17,
      l10wt > 0.5   ~ 16,
      l10wt > 0     ~ 15,
      l10wt > -0.5  ~ 14,
      l10wt > -1    ~ 13,
      l10wt > -1.5  ~ 12,
      l10wt > -2    ~ 11,
      l10wt > -2.5  ~ 10,
      l10wt > -3    ~ 9,
      l10wt > -3.5  ~ 8,
      l10wt > -4    ~ 7,
      l10wt > -4.5  ~ 6,
      l10wt > -5    ~ 5,
      l10wt > -5.5  ~ 4,
      l10wt > -6    ~ 3,
      l10wt > -6.5  ~ 2,
      l10wt > -7    ~ 1,
      l10wt > -7.5  ~ 0),
    # Weight ranges for log10() bins
    wt10_range = case_when(
      wt10bin == 21 ~ 10 ^ 3.5  - 10 ^ 3,
      wt10bin == 20 ~ 10 ^ 3    - 10 ^ 2.5,
      wt10bin == 19 ~ 10 ^ 2.5  - 10 ^ 2,
      wt10bin == 18 ~ 10 ^ 2    - 10 ^ 1.5,
      wt10bin == 17 ~ 10 ^ 1.5  - 10 ^ 1,
      wt10bin == 16 ~ 10 ^ 1    - 10 ^ 0.5,
      wt10bin == 15 ~ 10 ^ 0.5  - 10 ^ 0,
      wt10bin == 14 ~ 10 ^ 0    - 10 ^ -0.5,
      wt10bin == 13 ~ 10 ^ -0.5 - 10 ^ -1,
      wt10bin == 12 ~ 10 ^ -1   - 10 ^ -1.5,
      wt10bin == 11 ~ 10 ^ -1.5 - 10 ^ -2,
      wt10bin == 10 ~ 10 ^ -2   - 10 ^ -2.5,
      wt10bin ==  9 ~ 10 ^ -2.5 - 10 ^ -3,
      wt10bin ==  8 ~ 10 ^ -3   - 10 ^ -3.5,
      wt10bin ==  7 ~ 10 ^ -3.5 - 10 ^ -4,
      wt10bin ==  6 ~ 10 ^ -4   - 10 ^ -4.5,
      wt10bin ==  5 ~ 10 ^ -4.5 - 10 ^ -5,
      wt10bin ==  4 ~ 10 ^ -5   - 10 ^ -5.5,
      wt10bin ==  3 ~ 10 ^ -5.5 - 10 ^ -6,
      wt10bin ==  2 ~ 10 ^ -6   - 10 ^ -6.5,
      wt10bin ==  1 ~ 10 ^ -6.5 - 10 ^ -7,
      wt10bin ==  0 ~ 10 ^ -7   - 10 ^ -7.5
    ))  



# Summarise bins to get counts and normalizeed counts
bin_summs <- nefsc_stratified %>% 
  group_by(svvessel, lenbin2, wtbin, wt10bin, len_range, wt_range, wt10_range) %>% 
  summarise(n_group = sum(numlen_adj),
            .groups = "keep") %>% 
  ungroup() %>% 
  mutate(l_norm     = n_group / len_range,
         wtbin_norm = n_group / wt_range,
         wt10_norm  = n_group / wt10_range,
         `min length (cm)`  = as.character(round(len_range, 2)),
         `min weight (kg)`  = as.character(round(wt_range, 2)),
         `log10 min weight` = as.character(round(wt10_range, 2)))


# Use factor level formatting to put meaningful labels on the bins
sort_numeric_chars <- function(num_as_text){
  as.character(sort(as.numeric(unique(num_as_text))))
}

len_levels  <- sort_numeric_chars(bin_summs$`min length (cm)`)
wt_levels   <- sort_numeric_chars(bin_summs$`min weight (kg)`)
wt10_levels <- sort_numeric_chars(bin_summs$`log10 min weight`)


# apply factor levels
bin_summs <- bin_summs %>% 
  mutate(
    `min length (cm)` = factor(`min length (cm)`, levels = len_levels),
    `min weight (kg)` = factor(`min weight (kg)`, levels = wt_levels),
    `log10 min weight` = factor(`log10 min weight`, levels = wt10_levels)
  )

Length Bins

# log() length bins
l1 <- ggplot(bin_summs, aes(x = `min length (cm)`, y = n_group, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x = "Minimum Length (cm)", y = "Abundance") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

l2 <- ggplot(bin_summs, aes(`min length (cm)`, l_norm, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x = "Minimum Length (cm)", 
         y = "Normalized Abundance", 
         caption = "Bin numbers increase at increments of 0.5 in log(length).") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

l1 | l2

log(weight)

# log() weight bins
wt1 <- ggplot(bin_summs, aes(`min weight (kg)`, n_group, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x = "Minimum Weight (kg)", y = "Abundance") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

wt2 <- ggplot(bin_summs, aes(`min weight (kg)`, wtbin_norm, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x =  "Minimum Weight (kg)", 
         y = "Normalized Abundance", 
         caption = "Bin numbers increase at increments of 2 in log(weight).") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))


wt1 | wt2

log10(weight)

#l10 weight bins
wt10_1 <- ggplot(bin_summs, aes(`log10 min weight`, n_group, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x = "Minimum Weight (kg)", y = "Abundance") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))


# abundance normalized by bin width
wt10_2 <- ggplot(bin_summs, aes(`log10 min weight`, wt10_norm, fill = svvessel)) +
    geom_col(show.legend = F) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~svvessel) +
    scale_fill_gmri() +
    labs(x = "Minimum Weight (kg)", 
         y = "Normalized Abundance", 
         caption = "Bin numbers increase at increments of 0.5 in log(length).") +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5))

wt10_1 | wt10_2

Annual

ann_summs <- nefsc_stratified %>% 
  group_by(est_year, svvessel, lenbin2, wtbin, wt10bin, len_range, wt_range, wt10_range) %>% 
  summarise(n_group = sum(numlen_adj),
            .groups = "keep") %>% 
  ungroup() %>% 
  mutate(l_norm = n_group / len_range,
         wtbin_norm = n_group / wt_range,
         wt10_norm = n_group / wt10_range,
         `min length (cm)`  = as.character(round(len_range, 2)),
         `min weight (kg)`  = as.character(round(wt_range, 2)),
         `log10 min weight` = as.character(round(wt10_range, 2)),
         `min length (cm)` = factor(`min length (cm)`, levels = len_levels),
         `min weight (kg)` = factor(`min weight (kg)`, levels = wt_levels),
         `log10 min weight` = factor(`log10 min weight`, levels = wt10_levels))




# Plotting l10 weight bins

# functions to reduce code repetition
plot_l10bins <- function(x){
  x %>%
  ggplot(aes(`log10 min weight`, n_group, fill = svvessel)) +
        geom_col(show.legend = F) +
        scale_y_continuous(labels = scales::comma_format()) +
        facet_grid(est_year ~ svvessel, scales = "free") +
        scale_fill_gmri() +
        labs(x = "Min Weight (kg)", y = "Abundance") +
        theme(strip.text.y = element_text(angle = 0),
              axis.text.x = element_text(angle = 90, vjust = 0.5))
}

# Normalized abundances
plot_l10bins_norm <- function(x){
  x %>%
  ggplot(aes(`log10 min weight`, wt10_norm, fill = svvessel)) +
        geom_col(show.legend = F) +
        scale_y_continuous(labels = scales::comma_format()) +
        facet_grid(est_year ~ svvessel, scales = "free") +
        scale_fill_gmri() +
        labs(x = "Min Weight (kg)", y = "Normalized Abundance") +
        theme(strip.text.y = element_text(angle = 0),
              axis.text.x = element_text(angle = 90, vjust = 0.5))
}




# List to store them
ann_plots <- list()


# Albatross
ann_plots$AL <- ann_summs %>% 
  filter(svvessel == "AL",
         est_year >= 2000) %>% 
  plot_l10bins()


# Bigelow
ann_plots$HB <- ann_summs %>% 
  filter(svvessel == "HB") %>%
  plot_l10bins()


# Albatross Normalized
ann_plots$AL_norm <- ann_summs %>% 
  filter(svvessel == "AL",
         est_year >= 2000) %>% 
  plot_l10bins_norm()

# Bigelow Normalized
ann_plots$HB_norm <- ann_summs %>% 
  filter(svvessel == "HB") %>% 
  plot_l10bins_norm()



# Assemble them
ann_plots$AL | ann_plots$HB

ann_plots$AL_norm | ann_plots$HB_norm

Spatial Patterns

For large regions like Georges Bank and the Gulf of Maine, what kind of patterns are we seeing.

# Load the strata
survey_strata <- read_sf(str_c(res_path, "Shapefiles/BottomTrawlStrata/BTS_Strata.shp"))  %>% 
  clean_names() %>% 
  filter(strata >= 01010 ,
         strata <= 01760,
         strata != 1310,
         strata != 1320,
         strata != 1330,
         strata != 1350,
         strata != 1410,
         strata != 1420,
         strata != 1490) 


# Key to which strata = which regions
strata_key <- list(
  "Georges Bank"          = as.character(13:23),
  "Gulf of Maine"         = as.character(24:40),
  "Southern New England"  = str_pad(as.character(1:12), width = 2, pad = "0", side = "left"),
  "Mid-Atlantic Bight"    = as.character(61:76))


# Assign Areas
survey_strata <- survey_strata %>% 
  mutate(
    strata = str_pad(strata, width = 5, pad = "0", side = "left"),
    strata_num = str_sub(strata, 3, 4),
    area = case_when(
      strata_num %in% strata_key$`Georges Bank` ~ "Georges Bank",
      strata_num %in% strata_key$`Gulf of Maine` ~ "Gulf of Maine",
      strata_num %in% strata_key$`Southern New England` ~ "Southern New England",
      strata_num %in% strata_key$`Mid-Atlantic Bight` ~ "Mid-Atlantic Bight",
    TRUE ~ "Outside Major Study Areas"
  )) %>% 
  select(finstr_id, strata, strata_num, area, a2, str2, set, stratuma, str3, geometry)


# Make trawl data an sf dataset
trawl_sf <- nefsc_stratified %>% st_as_sf(coords = c("decdeg_beglon", "decdeg_beglat"), crs = 4326)

Map of Trawl Regions

# Plot to check
ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(data = survey_strata, aes(fill = area)) +
  coord_sf(xlim = c(-77, -65.5), ylim = c(34, 45.75), expand = FALSE) +
  guides(fill = guide_legend(nrow = 2)) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank())

Map of Ecological Production Units

epu_sf <- ecodata::epu_sf

ggplot() +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  geom_sf(data = epu_sf, aes(fill = EPU)) +
  coord_sf(xlim = c(-77, -65.5), ylim = c(34, 45.75), expand = FALSE) +
  guides(fill = guide_legend(nrow = 2)) +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank())

Regional Summaries

# Just Area, all seasons
area_summs <- nefsc_stratified %>% 
  group_by(survey_area) %>% 
  summarise(
    season = "Spring + Fall",
    lw_biomass_kg = sum(sum_weight_kg, na.rm = T),
    n_stations = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length = weighted.mean(length, weights = numlen_adj),
    .groups = "keep"
  ) %>% 
  ungroup()

# Area x Season
seas_area <- nefsc_stratified %>% 
  group_by(survey_area, season) %>% 
  summarise(
    lw_biomass_kg = sum(sum_weight_kg, na.rm = T),
    n_stations = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length = weighted.mean(length, weights = numlen_adj),
    .groups = "keep"
  ) %>% 
  ungroup()


# Combine those two
summs_combined <- bind_rows(area_summs, seas_area) %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))

summs_combined %>% 
  mutate_if(is.numeric,round, 2) %>% 
  arrange(survey_area,season) %>% 
  knitr::kable()

survey_area	season	lw_biomass_kg	n_stations	lw_biomass_per_station	mean_ind_bodymass	mean_ind_length
GB	Spring	438226.3	2996	146.27	0.88	40.33
GB	Fall	450348.8	3124	144.16	0.73	37.69
GB	Spring + Fall	888575.1	6120	145.19	0.80	38.94
GoM	Spring	348424.0	3631	95.96	0.63	35.06
GoM	Fall	619849.5	3738	165.82	0.70	36.53
GoM	Spring + Fall	968273.5	7369	131.40	0.67	35.84
MAB	Spring	752939.2	2544	295.97	0.83	44.44
MAB	Fall	104074.9	2476	42.03	0.66	25.95
MAB	Spring + Fall	857014.1	5020	170.72	0.78	38.68
SNE	Spring	460508.7	2869	160.51	0.54	36.58
SNE	Fall	210873.7	2743	76.88	0.45	32.75
SNE	Spring + Fall	671382.4	5612	119.63	0.51	35.08

# Year x Area
area_summs_y <- nefsc_stratified %>% 
  group_by(est_year, survey_area) %>% 
  summarise(
    season                 = "Spring + Fall",
    lw_biomass_kg          = sum(sum_weight_kg, na.rm = T),
    n_stations             = n_distinct(id),
    lw_biomass_per_station = lw_biomass_kg / n_stations,
    mean_ind_bodymass      = weighted.mean(ind_weight_kg, weights = numlen_adj),
    mean_ind_length        = weighted.mean(length, weights = numlen_adj),
    expanded_abund         = sum(expanded_abund_s),
    expanded_lwbio         = sum(expanded_lwbio_s),
    expanded_biom          = sum(expanded_biom_s),
    .groups = "keep"
  ) %>% 
  ungroup()

Total Biomass

Total biomass in the figures below is derived from weight at length relationships and reflects the mean expected biomass of the fish caught using the observed distribution and frequency of lengths.

area_summs_y %>% 
  ggplot(aes(est_year, lw_biomass_kg)) +
    geom_line() +
    facet_wrap(~survey_area, ncol = 2) +
    scale_y_continuous(labels = scales::comma_format()) +
    labs(x = "", y = "Total Biomass (kg)")

CPUE

CPUE displayed below is the mean total biomass, derived from weight at length relationships, per station for each region. These values have not been weighted to reflect the difference in area between regions.

area_summs_y %>% 
  ggplot(aes(est_year, lw_biomass_per_station)) +
    geom_line() +
    facet_wrap(~survey_area, ncol = 2) +
    labs(x = "", y = "Adjusted Biomass per Station (kg)")

Area Expanded Abundance

Area expanded abundance is calculated by taking the catch per station rates for each length class of each species and extending that density to the total areas of the region. Catchability is assumed to be 1 for this calculation though it is certainly lower. This extension of catch rates is highly sensitive to large swings due to high variances and zero inflation of abundances.

area_summs_y %>% 
  ggplot(aes(est_year, expanded_abund/1000000)) +
    geom_line() +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~survey_area, ncol = 2) +
    labs(x = "", y = "Area Expanded Abundance (millions)")

Area Expanded Biomass

The expected biomass estimates for the survdat$biomass data and the expected biomass from L-W regressions differ wildly once the survey transitions over to the Henry Bigelow. Not sure at this point if that is because the coefficients are off, or if there are some specific stations that need to be investigated, but its quite a large difference:

fscs <- area_summs_y %>% 
  ggplot() +
    geom_line(aes(est_year, expanded_biom  /1000000, color = "Shipboard Weights")) +
    scale_y_continuous(labels = scales::comma_format()) +
    facet_wrap(~survey_area, ncol = 1) +
    scale_color_gmri(reverse = T) +
    labs(x = "", y = "Area Expanded Biomass FSCS (million kg)") +
    theme(legend.position = "bottom", legend.title = element_blank())

lw <- area_summs_y %>% 
  ggplot() +
    geom_line(aes(est_year, expanded_lwbio /1000000, color = "LW Regression Weights")) +
    scale_y_continuous(labels = scales::comma_format()) +
    scale_color_gmri(reverse = F) +
    facet_wrap(~survey_area, ncol = 1) +
    labs(x = "", y = "Area Expanded Biomass LW (million kg)") +
    theme(legend.position = "bottom", legend.title = element_blank())

fscs | lw

 

---

A work by Adam A. Kemberling

Akemberling@gmri.org